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Here we report a precise computer simulation study of the static critical properties of the two- 
dimensional g-states Potts model using very accurate data obtained from a modified Wang-Landau 
(WL) scheme proposed by Caparica and Cunha-Netto [Phys. Rev. E 85, 046702 (2012)]. This 
algorithm is an extension of the conventional WL sampling, but the authors changed the criterion 
to update the density of states during the random walk and established a new procedure to windup 
the simulation run. These few changes have allowed a more precise microcanonical averaging which 
is essential to a reliable finite-size scaling analysis. In this work we used this new technique to 
determine the static critical exponents j3, 7 , and iz, in an unambiguous fashion. The static critical 
exponents were determined as /3 = 0.10811(77), 7 = 1.4459(31), and u = 0.8197(17), for the g = 3 
case, and /3 = 0.0877(37), 7 = 1.3161(69), and u = 0.7076(10), for the g = 4 Potts model. A 
comparison of the present results with conjectured values and with those obtained from other well 
established approaches strengthens this new way of performing WL simulations. 


I. INTRODUCTION 

Monte Carlo (MC) simulations are ubiquitous in the 
field of statistical mechanics, especially for the study of 
phase transitions and critical phenomena fli- Since the 
historical work of Metropolis et al Q , the most outstand¬ 
ing task in this context is the pursuit of new and more ef¬ 
ficient algorithms to overcome long time scale problems. 
Since there are few problems in the field of interacting 
systems for which one can find an exact solution, MC 
simulations became an indispensable tool. This is due to 
the massively increasing in computational power and fur¬ 
ther due to the development of more efficient algorithms. 
More recently, such development focused on the extended 
ensemble method, where one uses an ensemble different 
from the ordinary canonical with a fixed temperature, as 
in the original Metropolis algorithm. To name a fill ex¬ 
amples we have the multicanonical method and the 
exchange Monte Carlo method (parallel tempering) Q. 
Particularly, during the last two decades, a multicanoni¬ 
cal MC algorithm known as Wang-Landau sampling 
has been at the forefront of interest Q and has proven 
to be a very powerful numerical procedure for the study 
of phase transitions and critical phenomena 

The original idea of the WL algorithm is to measure an 
a priori unknown density of states of a given system iter¬ 
atively by performing a random walk in energy space and 
sampling configurations with probability proportional to 
the reciprocal of the density of states, resulting in a “flat” 
histogram. Despite being a well-established numerical 
procedure, it is clear that some improvements on the 
algorithm are indeed necessary to overcome some lim¬ 
itations during the simulation run. The method itself 
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was subject to several studies and various improvements 
to it have been proposed [HHil- By its turn, the MC 
algorithm used in this work is an extension of the conven¬ 
tional WL where some few changes produce more reliable 
and precise results. 

Considering the aforementioned comments, the pur¬ 
pose of this paper is twofold. First, to present a numeri¬ 
cally simple and accurate procedure to halting a regular 
WL simulation run. This is accomplished with a method 
proposed in Refs. [Hill. Second, to apply this tech¬ 
nique to the square two-dimensional g-state Potts model 
and compute the static critical exponents for g = 3 and 
4 states, showing that this method is also a helpful tool 
to address the achievement of critical exponents, a pos¬ 
sibility barely explored in the literature, the exception 
being the important works of Malakis et al [2ll-l^|. In 
the following we will make use of a combination between 
finite-size scaling theory and cumulant methods to lo¬ 
cate and evaluate the extrema of various thermodynamic 
quantities and estimate the static critical exponents. 

The outline of this paper is as follows: In section II we 
define the model . In section III we define the simula¬ 
tion procedure. In section IV we describe the finite-size 
scaling analysis. The results are discussed in section V. 
Section VI is devoted to the summary and concluding 
remarks. 


II. g-STATES POTTS MODEL 

The Potts model, proposed by Potts in the early 
1950’s, has stood at the frontier of research in statisti¬ 
cal mechanics since its formulation. It is an extension of 
the two states Ising model to g > 2 states. In this model, 
to each lattice site is attached a spin variable ai (defined 
on each site i) which takes on integer values I,...,g. 
Adjacent sites have an attractive interaction energy — J 
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whenever they are equal or 0 otherwise. The Hamilto¬ 
nian of the q-states ferromagnetic model (J > 0) can be 
written as [1^ 

n = (1) 

<i,3> 

where S is the Kronecker <5—symbol, and the sum runs 
over all nearest neighbors of In the low temperature 
regime the system is ordered, becoming disordered as T 
increases. In 2D, for g < 4 the phase transition is of 
second-order and discontinuous if <7 > 5. A proper order 
parameter cj) is 

9= -1, 2) 

q-l 

where N^ax is the “volume” occupied by the spins of the 
state q of largest population and N = [l^ . 

III. ENTROPIC SIMULATIONS 

The Wang-Landau method @ is based on the fact 
that if one performs a random walk in energy space with 
a probability proportional to the reciprocal of the den¬ 
sity of states, a flat histogram is generated for the en¬ 
ergy distribution. Since the density of states produces 
huge numbers, instead of estimating g{E), the simula¬ 
tion is performed for S{E) = \ng{E). At the begin¬ 
ning we set S{E) = 0 for all energy levels. The random 
walk in the energy space runs through all energy levels 
from Emin to Emax with a probability p{E E') = 
min(exp [S'(A) — S{E')], 1), where E and E' are the en¬ 
ergies of the current and the new possible configurations, 
respectively. Whenever a configuration is accepted we 
update H{E') H{E') + 1 and S{E') S{E') + Fi, 

where Fi = In fi, fo = e = 2.71828... and f^+i = y/Jl 
{fi is the so-called modification factor). The flatness of 
the histogram is checked after a certain number of Monte 
Carlo steps (MCS) and usually the histogram is consid¬ 
ered flat if H{E) > 0.8(H), for all energies, where (H) 
is an average over energies. If the flatness condition is 
fulfilled we update the modification factor to a finer one 
and reset the histogram H{E) = 0. 

Recent works |17H^ have demonstrated that (a) in¬ 
stead of updating the density of states after every move, 
one ought to update it after each Monte Carlo sweep 
[2^ (this providence avoids taking into account highly 
correlated configurations when constructing the density 
of states); (b) WL sampling should be carried out only 
up to In/ = \nffinal defined by the canonical aver¬ 
ages during the simulations (this saves CPU time, dis¬ 
carding unnecessary long simulations); and (c) the mi- 
crocanonical averages should not be accumulated before 
In/ < Inf micro defined by a previous study of the mi- 
crocanonical averaging during the simulation (the ruled 
out WL levels in these averages correspond to a micro- 
canonical termalization, since the initial configurations 


do not match those of maximum entropy). The adop¬ 
tion of these easily implementable changes leads to more 
accurate results and saves computational time. They in¬ 
vestigated the behavior of the maxima of the specific heat 

C{T) = {{E-{E))^)/T^ (3) 

and the magnetic susceptibility 

x{T) = L^{{m-{m))^)/T, (4) 

where E is the energy of a given configuration and m 
is the corresponding magnetization per spin, during the 
WL sampling for the Ising model on a square lattice. 
They observed that a considerable part of the conven¬ 
tional Wang-Landau simulation is not very useful because 
the error saturates. They demonstrated in detail that in 
general no single simulation run converges to the true 
value, but to a particular value of a Gaussian distribu¬ 
tion of results around the correct value. The saturation 
of the error coincides with the convergence to this value. 
Continuing the simulations beyond this limit leads to ir¬ 
relevant variations in the canonical averages of all ther¬ 
modynamic variables. 

Zhou and Bhatt [l2| demonstrated that when / is close 
to 1 the relative error 5g/g = 6\ng scales as Vln /. Con¬ 
versely, in Ref. [l^ it was shown that this convergence 
indeed holds, but the final result falls in a Gaussian dis¬ 
tribution around the true value. In this work it is also 
noteworthy that the convergence described in Ref. [l^ 
for the 1/t entropic sampling, where the authors argue 
that the logarithm of the density of states converges as 
1/Vi, is not reflected in the canonical averages, since 
for long simulations different runs do not converge to a 
unique value, moreover the results take on an erratic be¬ 
havior. 

Ref. [1^ also proposes a criterion for halting the simu¬ 
lations. Applying WL sampling to a given model, begin¬ 
ning from /s, we calculate the temperature of the peak 
of the specific heat defined in Eq. ([3]) using the current 
g{E) and from this time forth this mean value is updated 
whenever the histogram is checked for flatness. When the 
histogram is considered fiat, we save the value of the tem¬ 
perature 7/(0) of the peak of the specific heat. We then 
update the modification factor /i+i = y/J) and reset the 
histogram H{E) = 0. During the simulations with this 
new modification factor we continue calculating the tem¬ 
perature of the peak of the specific heat Tc{t) whenever 
we check the histogram for flatness and we also calculate 
the checking parameter 

e=\Ta{t)-Ta{0)\. (5) 

If the number of MCS before verifying the histogram 
for flatness is chosen not too large, say 10,000, then dur¬ 
ing the simulations with the same modification factor 
the checking parameter e is calculated many times. If 
e remains less than 10“'* until the histogram meets the 
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flatness criterion for this WL level, then we save the den¬ 
sity of states and the microcanonical averages and stop 
the simulations. When one adopts this criterion for halt¬ 
ing the simulations, different runs stop at different final 
modification factors. 

Having at hand the density of states, one can calculate 
the canonical average of any thermodynamic variable X 
as 




( 6 ) 


where {X)e is the microcanonical average accumulated 
during the simulations and /3 = l/fc^r, where T is the 
absolute temperature measured in units of J/ks and ks 
is the Boltzman’s constant. 

In Ref. [3 it was also observed that two indepen¬ 
dent similar finite-size scaling procedures can lead to very 
different results for the critical temperature and expo¬ 
nents, which often do not agree within the error bars. 
The way to overcome this difficulty is to carry out 10 
independent sets of hnite-size scaling simulations. In 
the present work, for each of these sets and for each 
Potts model {q = 3 and q = 4), we performed simula¬ 
tions for L = 32, 36,40,44,48, 52,56,64, 72, and 80 with 
n = 24,24,20,20,20,16,16,16,12, and 12 independent 
runs for each size, respectively. The final resulting values 
for the critical exponents were obtained as an average 
over all sets. 


IV. FINITE-SIZE SCALING 

According to finite-size scaling theory (26l - [^ from the 
definition of the free energy one can obtain the zero field 
scaling expressions for the magnetization, susceptibility, 
and specihc heat, respectively, by 


(7) 


and it is expected that this expression is also exact for 
q = 3, although a rigorous proof of this assumption is 
still lacking | 1 ^ . 

Following Refs. Bmi we can define a set of thermo¬ 
dynamic quantities related to logarithmic derivatives of 
the magnetization: 


Vi = 4[m^] - 3[m'^], (12) 

V 2 = 2[to^] - (13) 

V 3 = 3[m^] - 2[m% (14) 

V 4 = (4[m] - [m'^])/3, (15) 

F 5 = (3H - [m3])/2, (16) 

Ve = 2[m] - [m^], (17) 

where 


[m”] = In 


9(m") 

dT 


Using Eq. © it is easy to show that 


(18) 


Vj^^\nL + Vj{tL^^‘') (19) 

for j = 1,2, ..., 6 . Since the critical temperature Tc is 
known for both models, at the critical temperature t = 0 
and the Vj are constants independent of the system size 
and we can estimate Xjv by the slopes of Vj calculated 
at Tc- And then, with the exponent v at hand, we can 
estimate the exponents /3 and 7 by the slopes of the log- 
log plots of Eqs. © and dH]) calculated at the critical 
temperature Tc- 


V. RESULTS 



( 8 ) 


( 9 ) 


where t = (Tc — T)lTc is the reduced temperature, and 
a, /3, and 7 are static critical exponents which should 
satisfy the scaling relation [ 2 ^ 

2 — a = dv = 215 + ^. ( 10 ) 

The critical temperature for the Potts model (for q > 
4) is exactly known as 


In all simulations we carried out, the microcanoni¬ 
cal averages were accumulated beginning from /y, we 
adopted the MCS for updating the density of states and 
the jobs were halted using the checking parameter e. In 
Fig. [ 1 ] we show the evolution of the temperature of the 
maximum of the specihc heat during the WL sampling 
beginning from /g for a single run with L = 52 and the 
evolution of logio(£) during the same simulation. One 
can see that at the last WL level the logarithm of e re¬ 
mains less than -4 indicating that the simulation can be 
stopped at the end of /15. 

According to Eq. CB) the critical adimensional tem¬ 
perature for the q = 3 Potts model is given by 


ksTc _ 1 

J ln(l -h yg) 


(11) = 0.994972861... (20) 

^ ’ J H 1 + V3) 
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FIG. 1. (color online). Upper panel: Evolution of the temper¬ 
ature of the maximum of the specific heat during the WL sam¬ 
pling, beginning from fg for a single run. The dots show where 
the modification factor was updated. Lower panel: Evolution 
of the logarithm of the checking parameter e during the same 
simulation. 


Evaluating the thermodynamic quantities Eqs. m- 
(Hzl) at this temperature and taking into account Eq. 
m, we are able to determine i by the slopes of the 
straight lines that we obtain with respect to InL. For 
each of these six slopes we calculate ly = l/(p) with 
Ah' = and take an average with unequal un¬ 

certainties ( 3 ^ over them. In Fig. [5] we present this set 
of lines. From the linear fits to these points we estimate 
that i = 1.20847(41), yielding v = 0.82759(98). Nev¬ 
ertheless these values represent the result of only one of 
the 10 sets of finite-size scaling simulations which were 
carried out. Initially we run over all sets calculating v in 
order to determine this exponent to the best precision. 

At this point we take a moment to discuss which pro¬ 
cedure should be adopted to calculate the mean value of 
these 10 results, a single averaging or an average with un¬ 
equal uncertainties. In order to investigate the behavior 
of the data under these two procedures, we grouped the 
five first sets in a large one and the last five in another 
large set. Taking the averages with unequal uncertainties 
we obtained v = 0.82272(26) and 0.81658(38), while if we 
take just single averages neglecting the error bars, we ob¬ 
tain V = 0.8230(19) and 0.8165(22), in each of these two 
large sets. One can see that the former procedure leads 
to unrealistic error bars, whereas the later yields results 
that intersect within ±2cr errors. We therefore adopt 
the single averaging here and in all the further calcula¬ 
tions. In Table U the fourth column displays the values 
obtained in each set and the final result in the last line: 
i/ = 0.8197(17). 

Next, with the critical exponent v accurately deter¬ 


mined, we can use Eqs. 0 -® to evaluate the expo¬ 
nents ^ and ^ by the slopes of the log-log plots. In 
Fig. [3] and Fig. |4] we show this finite-size scaling be¬ 
havior for each exponent, obtaining ^ = 0.1298(28) and 
^ = 1.753(16), respectively. We then calculate {3 = 
with Z\/3 = ^Av + vA^, and similarly for 7 and Aj, 
obtaining /3 = 0.1063(23), and 7 = 1.435(13). Again, 
these values were obtained at the first folder. In Ta¬ 
ble U we show the results for all sets and the best es¬ 
timates in the last line, yielding j3 = 0.10811(77), and 
7 = 1.4459(31). Finally using the scaling relation given 
by Eq. (uni we determined the exponent a = 2 — 2/3 — 7 
with Aa = 2Zi/3 -|- Z\ 7 . These results are also displayed 
in Table H] giving a = 0.3379(28). 



FIG. 2. (color online) Size dependence of Vj at the critical 
temperature. The slopes yield Xjv. 



FIG. 3. (color online) Log-log plot of size dependence of the 
magnetization at T, = 0.994972861. 


For the q = 4 Potts model the critical adimensional 
temperature is given by 

= 0.910239226... (21) 

J ln(l + Vi) 
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Q 

P 

7 

V 

Q 

P 

7 

V 


q = 3 Potts model 



q = 4 Potts model 


0.352(17) 

0.1063(23) 

1.435(13) 

0.82759(98) 

0.550(26) 

0.0836(71) 

1.283(12) 

0.7123(13) 

0.351(12) 

0.1063(31) 

1.4367(62) 

0.82364(52) 

0.541(18) 

0.0951(35) 

1.269(11) 

0.7045(13) 

0.340(14) 

0.1120(32) 

1.4361(83) 

0.82068(56) 

0.547(35) 

0.0855(74) 

1.282(20) 

0.70907(91) 

0.334(14) 

0.1044(37) 

1.4569(69) 

0.81688(60) 

0.539(24) 

0.0950(43) 

1.271(16) 

0.7085(12) 

0.350(16) 

0.1093(39) 

1.4315(91) 

0.82601(47) 

0.504(33) 

0.0907(74) 

1.315(19) 

0.70185(77) 

0.341(16) 

0.1099(23) 

1.440(12) 

0.81797(88) 

0.519(27) 

0.0895(50) 

1.302(17) 

0.7041(12) 

0.337(12) 

0.1099(30) 

1.4428(66) 

0.82410(79) 

0.538(24) 

0.0946(66) 

1.273(11) 

0.7085(13) 

0.336(12) 

0.1053(33) 

1.4535(62) 

0.81632(63) 

0.528(15) 

0.0898(28) 

1.292(10) 

0.70671(83) 

0.326(13) 

0.1094(27) 

1.4552(85) 

0.81162(52) 

0.521(23) 

0.0958(50) 

1.287(14) 

0.7059(14) 

0.329(10) 

0.1071(27) 

1.4563(50) 

0.81249(45) 

0.550(46) 

0.056(12) 

1.337(24) 

0.7074(12) 

0.3379(28) 

0.10811(77) 

1.4459(31) 

0.8197(17) 

0.5084(48) 

0.0877(37) 

1.3161(69) 

0.7076(10) 

TABLE 1. q 

= 3 and q = 4 Potts models: 

Ten finite-size scaling results for the exponents a, P, 7 , 

and u. The last line shows 

the average 

values over all the 

runs. 






Method 



Q 

P 

7 


V 

q = 3 Potts model 

Conjectured value [1^ 


1 

3 

1 

9 

13 

9 


5 

6 

Kadanoff variational RG [3^ 


0.326 

0.107 

1.460 

0.837 

Monte Carlo RG [3^ 


0.352 

0.101 

1.445 

0.824 

This work 



0.3379(28) 

0.10811(77) 

1.4459(31) 

0.8197(17) 

q = 4 Potts model 

Conjectured value [1^ 


2 

3 

1 

12 

7 

6 


2 

3 

Kadanoff variational RG [3^ 


0.488 

0.091 

1.330 

0.756 

Duality invariant RG [3^ 


0.4870 

- 

- 


0.7565 

This work 



0.5084(48) 

0.0877(37) 

1.3161(69) 

0.7076(10) 


TABLE II. Estimates of a, {3, 7 , and compared to results obtained with other techniques and conjectured values. 



FIG. 4. (color online) Log-log plot of size dependence of the 
susceptibility at Tc = 0.994972861. 


All the plots and finite-size scaling procedures are com¬ 
pletely analogous to those we described above for the 


q = 3 case. In Table U we display the results for the 10 
folders and our final estimates yielding a = 0.5084(48), 
P = 0.0877(37), 7 = 1.3161(69), and iz = 0.7076(10). 

Such large repetitious handling of data for obtaining all 
these canonical averages and finite-size scaling extrapola¬ 
tions were possible only by using shell scripting [36l - l40l |. 
This is an exceptional tool for those who work with sim¬ 
ulations. 

As a final discussion, we compare in Table |lT] our fi¬ 
nal estimates of the critical exponents to other well- 
established values. It is possible to see a good agreement 
for the q = 3 case, especially between those obtained by 
numerical means. For the q = A case, our results are be¬ 
low those of Refs. Notwithstanding our results 

for P, 7 and v are closer to the conjectured ones when 
compared to these approaches. This is a clear indication 
that our procedure of carefully handling very accurate 
data obtained by an entropic sampling simulation is a 
powerful and reliable technique. 
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VI. CONCLUSIONS 


In this work we explored the static critical behavior of 
q = 3 and g = 4 Potts models within a high-precision and 
refined Wang-Landau procedure. All results are in very 
good agreement with those obtained from other well es¬ 
tablished approaches. The most striking conclusion from 
our analysis, in our opinion, is that it is possible to obtain 
reliable and very precise calculations of critical exponents 
from WL sampling provided that the appropriate imple¬ 
mentations adopted in this work are made. Most impor¬ 


tant, the implementation of the present method remains 
as simple as the original idea of the WL algorithm. A 
further critical test of our algorithm would be provided 
by an analysis of the critical behavior of multi-parametric 
spin systems, which is a hard task for any conventional 
WL approach. 
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